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ABSTRACT 

We study the effects of a large-scale, ordered magnetic field in protoplanetary disks on Type I planet 
migration using a combination of numerical simulations in 2D and 3D and a linear perturbation 
analysis. Steady-state models of such disks require the inclusion of magnetic diffusivity. To make 
progress using ideal MHD, we focus on simplified field configurations, involving purely vertical {B^) 
and azimuthal (B^p) field components and a combination of the two. For each of the models we 
calculate the locations of the relevant resonances and of the turning points, which delineate the 
propagation regions of the MHD waves that transport angular momentum from the planet to the 
disk. We use both numerical and semianalytic methods to evaluate the cumulative back torque acting 
on the planet, and explore the effect of spatial gradients in the disk’s physical variables on the results. 
We conclude that, under realistic (3D) circumstances, a large-scale magnetic field can slow down the 
inward migration that characterizes the underlying unmagnetized disk — by up to a factor of ^ 2 
when the magnetic pressure approaches the thermal pressure — but it cannot reverse it. A previous 
inference that a pure-H^ field whose amplitude decreases fast enough with radius leads to outward 
migration applies only in 2D. In fact, we find that, in 3D, a pure-H^ disk undergoes a rapid transition 
to turbulence on account of a magnetorotational instability that is triggered by the planet-induced 
appearance of a weak B^ component. 

Subject headings: accretion, accretion disks - MHD - planet-disk interactions - protoplanetary disks 


1. INTRODUCTION 

Our current understanding of planets is that they form in 
the circumstellar disks from which their host stars origi¬ 
nally accreted most of their mass, and that, particularly 
in the case of giant, largely gaseous, planets, this hap¬ 
pens on timescales that are shorter than the dispersal 
times (typically a few million years) of these disks. The 
unanticipated finding of “hot Jupiters,” giant planets or¬ 
biting within ^10 stellar radii of their host stars, has 
indicated that such planets can move inward from their 
formation sites farther out in the protoplanetary disk. 
This outcome has been attributed to two possible mech¬ 
anisms. The first involves a highly eccentric motion, in¬ 
duced by gravitational interaction with another planet 
(or planets) or with a companion star, and eventual 
circularization through tidal interaction with the cen¬ 
tral star (e.g., iRasio'fc FordlflO^ iWu fc Murravl [200^ 
IWu fc Lithwickll2011ll . This process can take place after 
the disk has already dispersed. In contrast, the second 
mechanism, which envisions a slow inward drift while the 
planet maintains a quasi-circular orbit, requires the pres¬ 
ence of a gaseous disk. In the latter picture, the planet 
exerts a gravitational perturbation on the surrounding 
gas, and the back reaction of this perturbation on the 
planet can generate a net torque that causes the planet 
to m i grate (see, e.g..lLubow fc Ida] l20RT|jMey_fcNelson 

l20I2LlBaruteau fc Massetll2013l and iBaruteau et al.ll20I,‘l 

for recent reviews). When a planet’s mass is sufficiently 

^ present address: Astronomy Department, Adler Planetar¬ 
ium, Chicago, IL 60605, USA; auribe@oddjob.uchicago.edu 


small that one can neglect its effect on the disk structure, 
and if, in addition, its radial motion is slower than the 
effective viscous diffusion in the disk so that it can also 
be neglected, the planet is said to undergo Type I migra¬ 
tion. As was originally shown by iGoldreich fc Tremaind 
dnii), the response of the gas in this case remains lin¬ 
ear and can be represented by the effects of resonances: 
Lindblad resonances (LRs) that occur both inside and 
outside the corotation radius Tc (where the planet’s an¬ 
gular velocity Dp equals the angular velocity D of the 
gas), and a corotation resonance at Tc. It turns out 
that, when the disk is isothermal, the LRs dominate the 
torque. At the LR locations, the planet exchanges an¬ 
gular momentum with the disk through the launching of 
rotation-modified sound waves that combine to produce 
a spiral wake. The waves from the outer region are ex¬ 
cited in gas that rotates more slowly than the planet: 
they thus remove angular momentum from the planet 
and induce inward motion. The waves from the inner 
region have the opposite effect. In a Keplerian disk, the 
outer LR associated with a given azimuthal mode num¬ 
ber lies closer to the planet than the corresponding inner 
LR. Therefore, if the density varies smoothly with ra¬ 
dius, the interaction of the planet with the outer disk 
is stronger and the planet migrates inward. The drift 
timescale for Type I migration is inversely proportional 
to the planet’s mass Mp and is estimated to be shorter 
than the disk’s lifetime even for an Earth-mass planet at 
a distance of ^ I AU from a solar-mass star. The migra¬ 
tion timescale stops decreasing with increasing Mp when 
the planet becomes massive enough to clear a gap in the 
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disk through a nonlinear interaction, which implies that, 
so long as the local disk’s mass remains > Mp, giant 
planets move with the radial velocity of the surround¬ 
ing gas (Type II migration). However, Type I migration 
remains relevant to giant planet evolution in the con¬ 
text of the core-accretion model — curre ntly the favored 
scenario for gi ant-planet formation ('e.g.. IArmitagel[200^ 
iRafikovI 1200911 — in which such planets form through 
gas accretion onto comparatively small rocky cores after 
they exceed a threshold of a few Mq. Within this mod¬ 
eling framework, the predicted rapid inward migration 
of the rocky embryos makes it difficult to account for 
the observed spa tial distribution of giant planets (e.g., 
Ilda fc Linl[200^ . 

The sign and magnitude of the net torque acting on a 
planet are determined by the difference between the op¬ 
posing contributions from the regions inside and outside 
the planet’s orbit as well as by the corotation torque 
from the co-orbital region, and can depend sensitively 
on the physical properties of the disk in the vicinity of 
the planet. The aforementioned “rapid inward migra¬ 
tion” conundrum has motivated an exploration of more 
realistic disk models than the smooth, isothermal con- 
hguration that was originally studied. These investi¬ 
gations are still ongoing, although they have already 
yielded several promising clues. One notable result has 
been the finding that the corotation torque can become 
dominant when the disk is nonisothermal, and that the 
planet may, in fact, migrate away from the star if the 
disk has an outward-decreasing specific entropy and a 
sufficiently robust angular momentum transport mech¬ 
anism to keep this torque unsaturated (see the above- 
cited reviews). Another intriguing possibility that has 
been discussed in the literature is the occurrence of disk 
inhomogeneities where the predicted direction of planet 
migration is reversed (e.g.,lMasset et al.ll200^ Ilda &: LinI 
l2008bl : iHasegawa &: Pudritdl2011h . In addition, several 
studies have considered magnetic effects. Magnetic fields 
are expected to play a number of key roles in protoplan¬ 
etary disks. They may contribute significantly to the 
angular momentum transport that enables gas to reach 
the central star — either through magnetohydrodynamic 
(MHD) turbulence (likely induced by the magnetorota- 
tional instability [MRI]), which involves a disordered, 
small-scale field that transports angular momentum in 
the radial direction, or by means of an ordered, large- 
scale field, which can give rise to vertical transport (likely 
involving a centrif ugally driven win d) through the disk 
surfy es (see, e.g.. iBalbusI 120111 and lKonigl &: SalmeronI 
120111 for recent reviews). In this connection, it has been 
proposed that the radial motion of small (Mp < M^) 
planets on scales of a few AUs could be dominated by 


stochastic migration induced bv MRI turbulence (e.g.. 

Nelsonll2005 

; IJohnson et al.ll2006l lAdams & Blochll2009l: 

Uribe et al.l 

20111). Another relevant effect is the trunca- 


tion of the disk by a sufficiently strong stellar magnetic 
field, which could help stall inw ard-migrating p lanets be¬ 
fore they reach the central star (jLin et al .11 19961 1. Finally, 
a large-scale, ordered magnetic field that permeates the 
disk and is well coupled to the gas can have a direct ef¬ 
fect on planet migration by modifying the disk response 
to the planet’s gravitational potential. This change is 
due to the fact, in the presence of such a field, per¬ 


turbations are no longer associated with hydrodynamic 
(HD) sound waves but rather with either one of the three 
types of MHD waves: slow-magnetosonic (SMS), fast- 
magnetosonic (FMS), and Alfven. 

The direct magnetic effect on planet m igration was firs t 
considered in the semianalytic work of lTerau^ (1200311 . 
who carried out a linear perturbation analysis for the 
case of an isothermal, 2D disk (i.e., one in which there 
are no variations along the rotation axis) that is perme¬ 
ated by a purely azimuthal magnetic field . The results of 
that s tudy were confirmed numerically bv iFromang et aj] 
(|2005ll . It was found that the perturbations induced by a 
planet embedded in such a disk give rise both to modified 
LRs, which occur at the locations where the azimuthal 
velocity of the gas relative to the planet equals the phase 
velocity of FMS waves in the azimuthal direction, and to 
so-called magnetic resonances, which occur at the loca¬ 
tions where the relative azimuthal velocity matches the 
phase velocity of a SMS wave propagating along the field 
lineO While the excitation at the modified LRs again 
gives rise to the launching of density waves that prop¬ 
agate away toward either the inner or outer disk, the 
waves excited at a magnetic resonance propagate in a re¬ 
stricted region around it. Since the magnetic resonances 
are closer to the planet than the LRs, their contribution 
to the torque dominates over that of the Lindblad torque 
if the magnetic field is large enough. It was also found 
that there is no corotation resonance in this case but that 
nevertheless the whole region around corotation (and not 
just narrow zones around the magnetic resonances) con¬ 
tributes to the total torque0 Similarly to the LRs, in 
this case the outer magnetic resonances induce inward 
migration whereas the inner ones have the opposite ef¬ 
fect. This suggests that, if the field amplitude decreases 
sufficiently rapidly with radius r, inward migration could 
be slowed down or even reversed. This was indeed veri¬ 
fied to be the situation for the 2D, purely a z imuth al field 
case, where, for example, iFromang et al.l l|2005ll found 
(for a disk of uniform temperature and density in which 
the local thermal-to-magnetic pressure ratio is 4) that 
inward migration is reduced (resp., outward migration 
occurs) when the field scales as r~^ (resp., r“^). 

The magnetic field configuration adopted in the afore¬ 
mentioned studies is not expected to arise naturally in 
protoplanetary disks. These disks likely form in the 
gravitational collapse of rotating molecular cloud cores 
that are threaded by a large-scale interstellar magnetic 
field. This picture is consistent with the hourglass mag- 

^ As we note in Section 12.31 the locations of the LRs actually 
correspond to turning points, rather than to genuine resonances, 
of the perturbation differential equation. 

® The ITeroueml II2003I ) 2D linearized perturbation analysis was 
carried out in the ideal-MHD limit. Under these conditions, it 
is conceptually not surprising that the corotation resonance was 
found not to be present in the magnetic case. This is because the 
corotation torque can be interpreted in terms of the drag exerted 
by particles that move on horseshoe orbits in the vicinity of the 
planet (('e.g.. rWardlll99ir ). In the ideal-MHD limit, the particles are 
“frozen” onto the magnetic field lines, and, for a purely azimuthal 
field, cannot cross the radius of the planet on a horseshoe orbit. It 
has, however, been found that, as the field weakens and the horse¬ 
shoe motions (which act to suppress the magnetic resonances) are 
reestablished, the corotation torque acquires a new magnetic com- 
pone nt that could poten tially cause the planet to migrate outward 
fe.g.. muilet et al.ll2013l l. 
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netic field morphology exhibited by some of these cores, 
which can be attributed to the pinching and dragging 
of the field lines by th e inflowing matter (see, e.g., 
iKonigl fc Salmer'onll201l[ l. An open field of this type is 
conducive to the driving of disk winds, and, in its sim¬ 
plest representation, has an even symmetry about the 
disk midplane (corresponding to the radial and azimuthal 
field components changing sign across the midplane even 
as the vertical component remains nonzero there). In 
this case, the gas located in the vicinity of a planet at 
the midplane of the disk would be threaded by a purely 
vertical field rather than by a purely azimuthal field as 
was assumed in the above studies. As a first s t ep tq - 
ward generalizing the earlier results, iMuto et all (|2008[1 
examined the effect of a purely vertical field on planet 
migration through analytic and numerical calculations. 
They showed that, in 2D, the only effect of the field is 
to change the thermal sound speed to the (higher) FMS 
wave speed, which leads to a reduction in the magnitude 
of the torque acting on the planet in comparison with 
the HD case. They also demonstrated that, in 3D, both 
SMS and Alfven waves are excited in the disk. However, 
their study was carried out within the framework of th e 
shearing-sheet approximation fe.g.. lNaravan et ^119871 1. 
which does not distinguish between the regions outside 
and inside the planet’s orbital radius and therefore does 
not allow a determination of the differential (net) torque 
acting on the planet. Furthermore, their numerical cal¬ 
culations in the 3D case were performed in the limit of a 
thermal-to-magnetic pressure ratio ^ 1, which is not rel¬ 
evant to protoplanetary disks (except perhaps near the 
boundary with the stellar magnetosphere). 

In this paper and its companion (jBans et al.ll20I4L here¬ 
after Paper H) we carry out a systematic study of the 
effects of a multi-component, ordered magnetic field on 
Type I disk migration using a linear perturbation anal¬ 
ysis and ideal-MHD numerical simulations. As noted 
above, protoplanetary disks are expected to have a sig¬ 
nificant radial field component on account of the inward 
gravitational pull of the central star. In the presence of 
the strong azimuthal shear that characterizes Keplerian 
disks, this component would be rapidly wound up into a 
strong azimuthal component that would get out of equi¬ 
librium unless the disk is diffusive. While real disks likely 
possess extended regions that have a sufficiently large 
diffusivity to justify a quasi-equilibrium treatment, the 
numerical modeling of such disks, especially if vertical 
angular momentum transport involving a centrifugally 
driven wind is taken into account, is a challenging and 
costly endeavor. We defer simulating planet migration in 
this general case to a future publication and concentrate 
in this paper on gaining physical insights into this process 
through a series of model problems that can be consid¬ 
ered within the framework of ideal MHD. Specifically, we 
investigate the cases of a purely vertical and a purely az¬ 
imuthal field, as well as the case of a combined vertical 
and azimuthal field, in 2 and 3 dimensions. Our general 
linearization formalism is described in Paper H, where it 
is illustrated with semianalytic results for a purely verti¬ 
cal field as well as for the case where an azimuthal field 
whose amplitude increases with height (from zero at the 
midplane) is present, alone or in conjunction with a ver¬ 
tical field (with the latter case enabling vertical angular 


momentum transport and thereby mimicking a real wind¬ 
driving disk). In Section[^of this paper we briefly review 
the linearization procedure and then present analytic and 
semianalytic results for the case of a field with vertical 
and azimuthal components that are uniform with height 
but variable with radius. Our numerical simulations are 
presented and discussed in conjunction with the results 
of the linear analysis in Section [3l with further discus¬ 
sion given in Section |4l We summarize our findings in 
Section [S] 

2. LINEAR ANALYSIS 

In this section we perform a linear perturbation anal¬ 
ysis to calculate the torque exerted by the disk on the 
planet. We expect the analytic results to complement 
the numerical simulations and to enable a fuller under¬ 
standing of the disk migration process in the presence of 
a large-scale, ordered magnetic field. For the vertical and 
azimuthal fields we investigate, we present the locations 
where the density perturbation in the disk is divergent 
(the resonances) and also the locations where the planet 
can excite waves in the disk (the turning points). A de¬ 
tailed general formulation for calculating the effect of a 
planet on a magnetized disk is presented in Paper 110 As 
we noted in S ection [H a nalyt ic studies were previously 
performed by iTerauemI (|2003[ ) for the case of a p urely 
az imuthal field in 2D (see also lFromang et ffll2005h and 
bv iMuto et al.l (1200811 for a purely vertical field in 2D and 
3D. 


2.1. Resonances 

We first review the basic formulation, concentrating on 
the case where the field has only vertical and azimuthal 
components (in an inertial nonrotating cylindrical coor¬ 
dinate system {r, </?, z}). In the ideal-MHD limit, the 
conservation equations for momentum and mass, and the 
induction equation, take the form 


dv 

- + (v.V)v 


=-VP-pVV' + —(VxB)xB , (1) 
47r 


^+V-(pv)=0, (2) 


—+Vx(vxB)=0, (3) 

where p, v, and P are the gas density, velocity, and 
pressure, respectively, ^|J is the gravitational potential, 
and B is the magnetic field. The equilibrium velocity 
is taken to be vq = (0,rDo,0), where D is the angu¬ 
lar velocity and the subscript ‘0’ denotes an equilibrium 
value. We assume that the disk is geometrically thin 
and rotates at nearly the Keplerian angular velocity; i.e., 
Do ~ Dk. The equilibrium magnetic field is given by 
Bq = (0, i? 2 pr“‘?*), where we henceforth nor¬ 

malize r by the radial location rp of the planet (subscript 

^ In Paper II we adopt the convention of previous semianalytic 
treatments and evaluate the torque that the planet exerts on the 
disk, which has the opposite sign to the torque exerted by the disk 
on the planet that we calculate in this paper. 























4 


Uribe, Bans, & Kbnigl 


‘p’). Unless noted otherwise, we take all equilibrium 
quantities to be uniform in z. The presence of the planet 
induces perturbations in these variables through its grav¬ 
itational potential, given by Equation (l20ll below. Work¬ 
ing in Fourier space in the t, tp, and z directions, we write 
all perturbed disk variables as 5X = 
with bj = mflp and m a nonnegative integer. The 
planet’s angular velocity idp is equal to the Keplerian 
angular velocity at r^. 

We dehne the Lagrangian displacement, in terms of 
the perturbed velocities as 

u' = ima^r , (4) 


v'^ = ima^^ , (6) 

where cr = rip — flp. Using B' = V x x B) 
(jChandrasekhar fc FermilflQb.lH and the above definitions 
for the linearized induction equation becomes 

B' = , (7) 

- ikzB^pr~‘>'^^^ - B^pr~'^'<’ ^ , (8) 

B' + {q^ - 

. (9) 

Linearizing Equations o - ®, plugging into Equa¬ 

tions O-®, and rearranging the resulting terms, yields 
a second-order differential equation for 

+ = Si{r)^ + So{r)'ilj' . 

( 10 ) 

When the leading-order coefficient A 2 vanishes, the dif¬ 
ferential equation becomes singular: the locations where 
this happens represent resonances. In the case of a field 
that has vertical and azimuthal components, the leading- 
order coefficient become^ 

A 2 = [m^a'^ - {vAzkz + vavK)'^] (H) 

X + c^) - c^{vAzkz + VA^k^f] 

X [{vAzkz + i^A^k^Y - m^cr^] 

-|-m^cr^(uA^^ - . 

^ The other coefficients in Equation unii cannot be written out 
in such a relatively sim ple form. All the functions that appear as 
coefficients in Equation (llOl l {Aq, Ai, A 2 , So, and ^i) are, however, 
provided as online supplementary material in the form of both a 
Mathematica notebook and a PDF document. 


Here c is the isothermal speed of sound (i.e., = P/p), 

and vaz, VAip are the Alfven speeds associated with the 
z and (fi components of the field, respectively (i.e., = 

B'/^/Anpo, v\^ = Bl^p/Anpo). For notational consistency 
we set k,p = m/r, and for the sake of brevity we use 
= '*^Az + ^A(^ “ k1 + k/p. The numerator of 

Equation CD vanishes in two different cases. For an 
arbitrary value of m, the first one applies when either 
kz ^ 0 or both vaz and VA,p are nonzero. It reads 

m^cr^ = {vAzkz + VA^k^p)"^ , (12) 

which indicates a resonance where the (Doppler shifted) 
forcing frequency ma equals the frequency of an Alfven 
wave. We therefore refer to this resonance as the Alfven 
resonance (AR). The second case is 


C^ivAzkz + VAvk^p)"^ 


'^Az ' 


■ 


(13) 


which indicates that there is another resonance where 
the forcing frequency matches that of a slow MHD wave; 
we refer to this second resonance as the magnetic reso¬ 
nance (MR). These four resonances (inner/outer AR and 
inner/outer MR) of a Bz, + By, fi e ld con figuration were 
previously discussed bv iFu &: Tail (|201ir) in the context 
of black-hole accretion disks. 


In the respective limits. Equation dm yields the pure- 
Bz and pure-R,^ cases that were considered in previous 
work. In particular, 


A2{Bp ^ 0) = (14) 

(mV^ - VAz^D/m'^C^'^ivAz + c^) " <^'>^Az^l) 
—c^k'^{v\^k1 — iAz^"^ ~ 

The numerator of Equation (03 vanishes at the Alfven 
and magnet ic resonance s that were previously found for 
this case bv iMuto et al.l (|2008[ 1. We note that these au¬ 
thors identified a third, “vertical” resonance from the 
vanishing of the denominator of Equation m in the 
limit kp ^ 0 (which corresponds to the WKB approxi¬ 
mation employed in their derivation). However, accord¬ 
ing to the definition given in Section 12.31 this is techni¬ 
cally a turning point and not a resonance. In the case of 
a pure-Ry, field in the 2D limit (kz = 0), Equation dill) 
gives, after factors of — v\^kp) are canceled from 

its numerator and denominator. 


A2{Bz = 0,kz = 0) = 


g’' (^Ay + C ) - C ^aAv 


(15) 

This resu l t repr oduces the MR identified for this case by 
iTerauenJ (|2003f ). In the limit kz = 0, Equations (IT^ and 
(fT3)l imply that the locations of the AR and the MR are 
independent of m, and, for a geometrically thin disk, are 
given (to order h) by 


r-AR = Tc ± 



2hp 

?'MR = "re ± , = , 

3\/l + Ppp + Pipp/ Pzp 


(16a) 

(16b) 
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where /3 is the ratio of the squares of the sound and 
Alfven speeds (i.e., j3z = ^^ = 

c/(O k?’p). 0 For a purely azi muthal field, Eq uation (Il6b|) 
reproduces equation (40) in iTerau^ (|2003fl . 


The locations of the re sonances ar e in general obtained 
with respect to ('e.sf.. lWai^ll997[ l. so we are interested 
in tracking this radius. Using the radial component of 
Equation o for a locally isothermal disk in which the 
midplane density and the sound speed scale with radius 
as po(z = 0) = Ppr““ and cq = Cpr~°'% we obtain 

1 = —- -({a + 2as)r^^°“‘+ (17) 

ri Tc V 


r 


a 

c 



(gy -2q^ 

/3^P ^ 


where the P parameters are evaluated at the midplane. 
The dependence of on is illustrated in Figure [U 
using the disk models employed in Figure |S] below. 



Fi g. 1. — Location of the normalized corotation radius (Equation 
(|17l> ') as a function of the azimuthal field-strength parameter |3^p in 
a self-similar disk model (a = 3/2, as = 1/2, = 5/4; 

note that the (3 parameters are spatially constant in this model) 
for a pure-Bt^ field oo; solid line) and for a Bz + B^p field 

configuration with I3z = (dashed line). The diamond symbols 
mark the values of rc for the two models in Figure [8] 

2.2. Resonances in a Strict-2D Limit 

In order to compare the results of the analytic formula¬ 
tion with those obtained from our 2D numerical simula¬ 
tions (in which no vertical motion is allowed) we derive 
the leading-order coefficient of Equation for the case 
where the condition ?;(, = 0 is enforced. We note that this 
procedure is not fully self-consistent when both vertical 
and azimuthal equilibrium field components are present 
since the perturbed magnetic force term in Equation m 
has a vertical component in this case that will induce 
vertical motion. We nevertheless proceed with this for¬ 
mulation since it is useful for interpreting the behavior 
of the 2D simulation results. 


Under the above assumption, the leading-order coeffi¬ 
cient of Equation (|T0l) becomes 


^ 2,20 = 




^Az 


+ C2)- 


2 2 2 2 
^Az^Acp - 


^Az ' 


(18) 


This leads to the following condition for the MR: 


2 2 

r cr = -TT ^— 


■viz) 




^Az 


(19) 


T his equation take s the exact same form as equation (39) 
in iTeroueml ()2003[ 1. which gives the MR condition for a 
purely azimuthal field, if the sound speed, c, is replaced 
by the effective sound speed, CeS = \/+ vi^. There¬ 
fore, in this strict-2D limit, the effect of a vertical field 
component is just to increase the sound speed, which 
pushes the resonances and the turning points (see Sec¬ 
tion [TS]) away from rc. This behavi or of a vertical fiel d 
in this limit was first pointed out bv IMuto et al.l (1200811 . 

Interestingly, the implications of Equation (fTOl) are the 
opposite of those inferred from the = 0 limit of Equa¬ 
tion (fT^ . In the latter case, the addition of a vertical 
field to an azimuthal field pushes the M Rs closer to the 
corotation radius (see Equation (I16bl) l. The turning- 
point locations also differ for these two cases (see Fig¬ 
ure These results highlight the fact that a strictly 2D 
analysis may not capture the true response of the system. 
In the next subsection we return to the general case and 
consider the physical locations of the various resonances 
and of their associated regions of wave propagation. 

2.3. Turning Points 

The locations where the solution to Equation (fTOll goes 
from wave-like to evanescent are known as the turning 
points. At these radii the planet interacts most strongly 
with the disk. The size and radial position of the re¬ 
gions of wave propagation are important predictors of 
the rate and direction of planet migration. In general, 
the closer the wave propagation regions are to the planet, 
the stronger will the disk-planet interaction be. To pre¬ 
vent confusion in the ensuing discussion, we remind the 
reader that, by this definition, the Lindblad “resonances” 
are a ctually turning points (see iGoldreich fc Tremaind 
[197^ . 

We find the turning points by changing variables so that 
the homogeneous equation corresponding to Equation 
cni) has only second- and zeroth-order derivatives, and 
then solve numerically for the roots of the coefficient in 
front of the zeroth - order term — this procedure is de¬ 
tailed in iTerauer^ (|2003f) and is also described in Pa¬ 
per II. In the examples that we exhibit, each field com¬ 
ponent has Pp = 1 (so that, in cases where both B^p 
and B^p are nonzero, the total field is characterized by 
Pp = 0.5)0 We henceforth also adopt the value = I 


^ In the limit of a thin disk, h (<C 1) is equal to the (normalized) 
tidal density scale height, which, for a vertically uniform magnetic 
field, is a good measure of the effective half-thickness of the disk. 


‘ Note, however, that the linear analysis results are based on 
a sound speed (cp = 0.03) that is slightly lower than the value 
(cp = 0.1) used in the numerical simulations. The reason for this 
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for the radial power-law exponent of the equilibrium az¬ 
imuthal field component, which corresponds to a “force 
free” configuration (no radial force associated with this 
field component; {l/r)B^pd{rB^p)/dr = 0). 

Figure [5] shows the turning points and regions of wave 
propagation for four different field configurations in the 
2D limit. In the case of a purely vertical field (panel (a)), 
there is only one set of turning points, corresponding 
to the effective Lindblad resonances@ At high values 
of m, the locations of the effective Lindblad resonances 
deviate from those of the nominal ones due to a pres- 
sure effect, leading to a sharp cutoff in the tor que (e.g., 
iGoldreich fc Tremainelll98^ lArtvmowiczlllQ^ . A finite 
Bqz component increases the total pressure (or, equiva¬ 
lently, the effective sound speed Ces = \J(? + 
centuating this effect and causing the positions of the 
high-TO turning points to shift farther away from the cor¬ 
responding nominal resonances and thus remain farther 
away from the corotation radius (see Figure II.3, where 
the label ‘IF henceforth r efers to Pap e r II). Therefore, 
as was previously noted bv iMuto et al.l (|2008D . the effect 
of a vertical field in 2D is to reduce the magnitude of 
the torque exerted on the planet. We label the modified 
(by the presence of a magnetic field) effective Lindblad 
resonances by 

In the case of the purely azimuthal field shown in 
panel (b) of Figured additional turning points appear on 
the two sides of the MR, corresponding to regions where 
SMS waves can propagate. We label the turning points 
interior and exterior to the MR by Rm- and Rm+, re¬ 
spectively. In 2D, Rm- exists everywhere only above 
som e critical value of m that depends on Pip and h/r 
Isee lTeraueml[200^ . The presence of additional regions 
of wave propagation, which are closer to the planet than 
the regions where the FMS waves associated with the ef¬ 
fective Lindblad resonances propagate, suggests that the 
torque for this field configuration is stronger than for the 
pure-i?z case. When the field strength is increased, all 
the turning points {R m-, Rm-i-, and Rc-i-) shift away 
from Tc (see figure 1 in iTerau^ I2003D . which tends to 
decrease the torque. However, the amplitude of the SMS 
waves in the wave propagation regions surrounding the 
magnetic resonances goes up when the field gets stronger, 
which has the opposite effect. It turns out that the lat¬ 
ter effect dominates, resulting in a larger torque for a 
stronger field in this case (in contrast with the behavior 
of a pure-i ?2 disk). 

As was shown in Section [521 when one enforces the con¬ 
straint v'^ =Q for a disk with both Bp and B^ field com¬ 
ponents, the magnetic resonance condition is the same 
as for a pure-R^ field except that c is replaced by the 
(larger) effective sound speed Ceff (see Equation (IT9)) L 
The same effect is found in the behavior of the turning 
points for this field configuration, shown in panel (c) of 

difference is that a comparatively small value of c lower s th e com¬ 
putational requirements on the integration of Equation ra by re¬ 
ducing the size of the integration domain, but that a larger sound 
speed facilitates resolving the locations of the resonances in the 
simulations. 

® The appellation “effective” is used to distinguish the turning 
points from the nominal Lindblad resonances, defined by m(f2 — 
Op) = ±K, where k is the disk’s epicyclic frequency. 


Figure [51 all the turning points are shifted away from 
the corotation radius, as expected on the basis of such a 
substitution. This suggests that the disk would behave 
similarly to the pure-R,^ case, but that the net torque 
would be reduced (resulting in slower inward migration). 
When the restriction u), = 0 is removed (panel (d)), one 
finds an additional set of turning points, with a corre¬ 
sponding new region of wave propagation. These turning 
points straddle the Alfven resonance position, which, like 
the MR position in this case, is independent of m (see 
Equation (fT6)l 'l. We label the interior and exterior Alfven 
turning points by Ra- and Ra+, respectively. We also 
find that the MR, Rm- , and Rm+ loci now lie closer to 
the corotation radius than in the = 0 case. This fact, 
along with the appearance of new regions of wave prop¬ 
agation, lead to a stronger torque for the u' 0 setup, 
although the net torque for this case is still weaker than 
for the pure-R,^ configuration (see left panel in Figure |7| 
below). 

As is clear from our results in Section [2Tl in 3D (fez 0) 
both the magnetic and the Alfven resonances, as well 
as their associated turning points and wave-propagation 
regions (SMS and Alfven, respectively), are present for 
either the pure-Rz, pure-R,^, or their combined field con¬ 
figuration. The latter two cases are shown in Figure O 
For the pure-R,^ disk shown in panel (a) of this figure, 
the innermost turning point {Rm-) is only present above 
a certain critical value of m, a behavior also seen in 2D 
(see panel (b) of Figure [2]). However, we find that in 3D 
this value also has a slight dependence on the vertical 
wavenumber, becoming larger as fez is increased. The 
fact that wave propagation in the immediate vicinity of 
the planet is possible for low values of m suggests that 
these azimuthal mode numbers would be the dominant 
contributors to the torque on the planet in this case. In 
contrast with this behavior, the innermost evanescent re¬ 
gion extends over all values of m for the Bp + Rz field 
configuration shown in panel (b) of Figure |3l This struc¬ 
ture is evidently induced by the vertical field component, 
since it is similar to that of a 3D pure-Rz field (see Fig¬ 
ure H.2). However, in the latter case the innermost turn¬ 
ing points disappear for su fficiently low valu es of 
(see Section H.2.5 as well as iMuto et al.ll200^ — we have 
confirmed that, when this product is small, the behavior 
of Rm- resembles that of the innermost turning points 
in panel (a) of Figure [3] 

3. NUMERICAL SIMULATIONS 
3.1. Numerical Setup 

Our numerical simulations are performed with version 
4.0 of the PLUTO code ([Mignone et al.l I2007D . We 
employ the HLLC and HLLD approximate Riemann 
solvers for, respectively, HD and MHD computations, 
and use a second-order picewise parabolic spatial in¬ 
terpolation method and a second-order Runge-Kuta 
solver for time integrations. The Constrained Transport 
method is adopted in the MHD simulations to preserve 
a divergence-free magnetic field. Our calculations are 
carried out in the ideal regime (no explicit viscosity or 
resistivity). We use polar coordinates r = {r,ip) in the 
2D simulations, and cylindrical coordinates r = (r, tp, z) 
in 3D. 




















Type I Planet migration in a Magnetized Disk. I. 


7 




Fig. 2.— Turning-point locations as a function of the azi mut hal wavenumber m in 2D for disks with the following field confi gura tions: 
(a) pure Bz^ (b) pure Bip, (c) Bz + with = 0 fSection l2.2ll . (d) Bz Bip in the 2D limit {kz = 0) of the 3D model (Section [23 . The 
radial distance is normalized by the planetary orbital radius rp, and each field component is characterized by /3p = 1. The shaded regions 
represent zones of wave evanescence, the dash-dotted lines mark the locations of the magnetic (MR) and Alfven (AR) resonances, and the 
solid lines show the locations of the turning points (whose labels are described in the text). 


The planet is assumed to move in a fixed, circular orbit, 
and is included in the simulation as an extra gravitational 
potential felt by the circumstellar disk. This potential 
has the form 


-ijj'ir) 


GMp 

(|r-rp|2 +£2)1/2 


( 20 ) 


where G is the gravitational constant and e is a soften¬ 
ing parameter (whose choice is discussed in ApoendixlXl). 
The planet’s mass is taken to be a fraction 1.5 x 10“® of 
the stellar mass M*, and, using the same normalization 
as in Section [21 the magnitude rp of its orbital radius 
is set equal to 1. We carry out the simulations in the 
frame of reference of the star, neglecting the (small) stel¬ 
lar motion around the center of mass of the star-planet 
system. For the 2D simulations, the planet is assumed 
move in the x — y plane, and the computational domain 
is given by r G [0.4, 2.5] and Lp G [0,27r]. In this case 


only the radial and azimuthal equations are solved, and 
all vertical gradients {d/dz) are taken to be zero. In 
2D, the number of cells in the computational domain 
is {Nr, Nip) = (512,1024) and the radial resolution is 
Ar = 0.004 (approximately 1/10 of the radial distance 
from the planet to the MRs). For the 3D calculations, 
the radial and azimuthal domains are the same, whereas 
the vertical domain is given by z G [—0.13,0.13]. In 
3D, the number of cells in the computational domain is 
{Nr, Nz, Nip) = (512,64,512) and the vertical resolution 
is equal to the radial one (Az = 0.004). We calculate the 
z component of the specific torque T exerted by the disk 
on the planet (hereafter the torque) by integrating over 
the volume V of the entire computational domain. 


Tz = G p{r) 


(i-p X r). 


(|r-rp|2 


2)3/2 


dV, 


( 21 ) 


where (rp x r)^ = (rp x r) • e^, and is the Cartesian 
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(b) 


Fig. 3.— Turning-point locations as a function of the azimuthal wavenumber m in 3D for disks with the following field configurations: 
(a) pure Btp, (b) B^p + Bz, with each field component being characterized by 0p = l.The vertical wavenumber is given by kzh = 1.56 in 
both panels, with h = c/Dp being the nominal tidal scale height. Notation is the same as in Figure[2] 


unit vector in the z direction. 

We refer to the running time average of the specihc 
torque as the cumulative torque. Torque is given in units 
of evaluated at rp. 

Initial Conditions 

The initial density prohle of the disk is given by po = 
Ppr~°‘. The adopted equation of state of the disk gas is 
P = pc^, where c is the isothermal sound speed (which 
we also take to have a power-law dependence on radius: 
c = Cpr“““). The tidal scale height is normalized by 

hp = 0 . 1 . 

Continuing to follow our analytic formulation, we take 
the initial field to be either purely azimuthal Bq = 
(0, 0), purely vertical Bq = {0,0, Bzpr~'^’‘), or 

with both components Bq = (0, , B^pr'?*). We 


again parametrize the midplane field amplitudes by the 
/3 parameters: B,pp = yJpPp/Pipp, B^p = yJpPp/Pzp, 
where p = 47r, and take /3,^p = /3^p = 1 as the default 
values. (The motivation for these choices is that, in a 
well-coupled, wind-driving disk, decreases with height 
from a value > 1 at the midplane, whereas (3ip, which 
formally diverges at the midplane for an even field sym- 
met ry, typically satisfies 0z < 10 at the disk surface; 
e.g., I'Wardle 

Units 

The density is giyen in physical units of Punit = 
10“^^ g cm“^, distances in units of runit = 1 AU, and 
yelocities in units of Uunit = \/GMQ/runit- Magnetic 
fields are giyen in units of Bunit = \/pPunitr^unif 

Boundary conditions 
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For the 2D simulations, the inner and outer radial bound¬ 
aries have zero-gradient boundary conditions for all phys¬ 
ical quantities. We impose < 0 at the inner radial 
boundary and n,. > 0 at the outer one. To prevent the de¬ 
velopment of high-frequency oscillations in the magnetic 
fie ld, we apply the wav e-damping procedure described 
in iFromang et al.l (|2005fl , which we implement by damp¬ 
ing radial and azimuthal velocity perturbations near the 
radial boundaries (r < 0.55 and r > 2.0) with an expo¬ 
nential function. The azimuthal boundary has a periodic 
boundary condition. 

For the 3D simulations, we use the same boundary con¬ 
ditions as in 2D in the radial direction, but we apply no 
damping to velocity perturbations during the initial ad¬ 
justment of the disk to a steady state. At the vertical 
boundaries, we implement the zero-gradient condition for 
all physical quantities whereas at the lower and upper z 
boundaries we impose Vz < 0 and Vz > 0, respectively. 

3.2. Pure-Bz Field in a 2D, Self-Similar Disk 

The first case we consider is that of a purely vertical field, 
which represents the midplane structure of a disk that is 
threaded by open magnetic field lines and possesses an 
even field symmetry with respect to the z = 0 plane (the 
generic model for a wind-driving d isk). This c a se wa s 
previously studied numerically by iMuto et all (|2008D . 
albeit under several restrictive simplifications (see Sec¬ 
tion [T]). For our simulations we adopt a = 3/2,as = l/2 
,and Qz = 5/4 — this choice of radial gradients is 
motivated by the scalings obtained from radially self- 
similar models of magnetized wind-d r iving disks ( e.g., 
iBlandford k. Pavn^ll982l : lK6nig]|[l989l : lTeitled[2011h . 

In the 2D limit, the magnetic term in the radial momen¬ 
tum equation {—(Bz/^i)dBz/dr) acts simply as an extra 
radial pressure force, so we expect the disk to behave as 
if it had a larger sound speed. As we discussed in Sec- 
tion[21 the effect of this on the effective LRs is to shift the 
locations of the turning points by a factor (1 -|- 
away from r^. This leads to a weaker coupling of the 
disk to the gravitational potential of the planet, which 
reduces the amplitude of the induced perturbations at 
these locations. Consequently, there is a reduction in 
the magnitude of the torque acting on the planet. In 
addition, for the adopted radial scalings, Equation Cl} 
shows that the corotation point is shifted inward (so that 
Tc < 1), the shift being larger for a stronger field (i.e., a 
lower /3z)- The effect of this inward shift is to increase the 
magnitude of the torque exerted by the outer disk while 
reducing the contribution (and potentially even changing 
the sign) of the torque exerted by the inner disk. The net 
torque on the planet is affected by both of these factors, 
the shift of the turning points and the shift of Cc, but in 
general the shift of the effective LRs away from will 
be the dominant effect (resulting in a lower torque from 
both the inner and outer disk regions). 

Our simulations confirm these predictions. Figured (first 
2 panels) shows the final (steady state) density structure 
in the r — ip plane for an HD and a pure-H^ disk. It 
demonstrates that the same basic spiral wake structure, 
arising from the constructive interference of waves ex¬ 
cited at the LRs, characterize both of these models. Fig¬ 


ure O in turn, shows the cumulative torques exerted on 
the planet as a function of time for the fiducial pure-52 
case {Pz = 1) as well as for disk models having fiz = 0.25, 
3, and oo (the HD limit). We present the separate contri¬ 
butions of the regions interior and exterior to the planet 
as well as their sum. The ability to evaluate the dif- 
fe rential torqu e disti nguishes our simulations from those 
of IMuto et al.l (|2008f) . which were carried out under the 
shearing-sheet approximation. We find that the planet 
migrates inward, as in both the HD case and the “pure 
Bz in a uniform disk” case studied semianalytically in 
Paper H (see Figure H.5). Our numerical results verify 
that the magnitude of the net torque decreases (implying 
slower inward migration) with increasing field strength. 
It is seen that the torque exerted by the inner disk is 
almost completely cut off for /3z < 1- In fact, for these 
low values of fdz it even contributes (albeit very weakly) 
with the same sign (< 0) as the outer disk. This seem¬ 
ingly paradoxical behavior is a reflection of the fact that 
the effect of the shifted LRs, which embody the disk’s 
response to the planet’s perturbing potential, is deter¬ 
mined by their location relative to Tc rather than Cp (see 
Section El]), and that for the \ow-j3z models the region 
between and rp becomes large enough to dominate 
the “inner” torque0 In a real wind-driving disk, how¬ 
ever, the magnetic pressure w ould not exceed the therma l 
pressure at the midplane le.g. lKonigl &: Salmeronl[20IIf) . 
For the lowest value of fiz that we simulate (0.25), the 
net torque is reduced to ^ 40% of the HD value. 

3.3. Bz + Btp Field Configuration in 2D 

We now generalize the pure-Hz disk model considered 
in the preceding subsection to the case where the field 
has both vertical and azimuthal components, which we 
take to be of equal magnitude. This corresponds to the 
field morphology near the surface of a wind-driving disk 
(although note that, in a real system, the vertical field 
component is typically stronger for a disk that is well- 
coupled throughout, and, in additio n, a radial field com¬ 
pone nt is also generally present; e.g.. lKonigl k. Salmeronl 
1201 If) . We still consider just the 2D limit. We present 
results for two models: a uniform disk (a = Ug = 0) 
with a “force-free” field configuration (qz = 0, q,p = 1), 
in which the contributions of both Bz and B,^ to the ra¬ 
dial magnetic force density vanish, and a self-similar disk 
(Qz = <}<p = 5/4, a = 3/2, Og = 1/2) akin to that con¬ 
sidered in Section [321 For reference, we also performed 
simulations of a pure-H^, disk model, complementing the 
analytic results for this field configuration given in Sec¬ 
tion [2] We do not discuss this case in de tail since s imilar 
calc ulations were pr e viousl y presented in lTeraue^ (j2003f ) 
and IFromang et ^ (|2005D . However, we show relevant 
plots for this model in Figures |S}iZl 

3.3.1. Force-Free Field in a Uniform Disk 

As we pointed out in Section 12.21 the numerical imple¬ 
mentation of the 2D limit is in general more constrain- 

^ As can be seen from Equation the shift in Tq depends not 
just on the value of 02 ^ but also on the radial power-law exponents 
of the density, sound speed, and magnetic field. The corotation 
radius is shifted to a smaller radius than rp when these quantities 
decrease with r. 
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Fig. 4.— Density in a uniform, 2D disk from simulations representing (from left to right) the HD (unmagnetized) case, a pure-B^ field, 
a pure-H(p field, and a Bz + B^p configuration. The planet is located at r = 1.0 and ip = tt. All models show the spiral wake induced 
by the Lindblad resonances, and disks with an azimuthal field component also exhibit magnetic resonances (manifested as the density 
perturbations that extend along the (p axis) both interior and exterior to the planet’s location. Each field component is characterized by 

0p = 1. 



Time (planet orbits) 


Fig. 5.— Cumulative torques on the planet in a 2D disk as a 
function of time for a pure-^z field configuration. The different 
lines, specified in the legend, present the contributions from the 
regions interior and exterior to the planet, as well as the net torque, 
for different values of the magnetic field-strength parameter /^z- For 
these runs, the values of the power-law exponents for the radial 
dependence of p, c, and Bz are a = 3/2, as = 1/2, and Qz = 5/4, 
respectively, corresponding to a self-similar disk. 

ing than the linear-analysis condition kz —>■ 0, and, for 
the field configuration we now consider, can be mimicked 
only by setting = 0 in the analytic derivation. Our 
discussion in Section 12.21 revealed that the presence of 
an azimuthal field component in a 2D disk introdnces 
magnetic resonances, and that the addition of a verti¬ 
cal field component shifts their locations away from rc 
(which, for the parameters adopted in this subsection, 
is eqnal to see Equation (fT7|) h The last 2 panels in 



Fig. 6.— Midplane density structure in the vicinity of the planet 
for two strict-2D field configurations. The dips in the plots mark 
the locations of the magnetic resonances. The parameters of the 
uniform, force-free disk model are a = 0, as = 0 , Qz = 0, q^p = 1, 
and each magnetic field component is characterized by = 1. 

Figure m demonstrate the appearance of the MRs when 
a component is present, and Figure [6] confirms the 
predicted outward shift in the resonance positions when 
a Bz component is added. In fact, according to Equation 
HU), for the chosen parameter values tmr = rp ± 0.047 
for the pure-R,^ case and tmr = rp±0.066 for the B^+B^ 
configuration. These compare well with the locations of 
the resonant features in the density plots shown in Fig¬ 
ure E The radial profile of the resonant features in each 
of these plots follows the structure of a tightly wound spi¬ 
ral perturbation that does not propagate radially, and so 
it changes with azimuth; in the example shown in Fig- 
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ure m the curves correspond to an azimuthal cut that is 
close to the planet’s location (0p — (j)cut = 0.09). 

Figure [7] shows the cumulative torque exerted on the 
planet, calculated semianalytically (left panel) after ob¬ 
taining the perturbed density from an integration of 
Equation (flOll (see Section II.3 for details), and numer¬ 
ically (right panel) from the simulation results. In the 
left panel, the torque is plotted as a function of the az¬ 
imuthal mode number m, whereas in the right panel it is 
shown as a function of time. Both derivations show that 
the net torque for either the pure-B,^ or the B^ + B^ field 
configuration is positive, corresponding to outward mi¬ 
gration. The result for the pure-B,^ with = 1 field in a 
uniform disk was previously obtained bv lTerauemI 1)200311 
(see top panel of figure 6 in that paper), and the fact that 
the same trend is exhibited by the B^ + B^ field configu¬ 
ration is not surprising in view of the turning-point anal¬ 
ysis in Section [7751 which indicated that the addition of 
a vertical field component weakens the torque but does 
not change the basic behavior of the system. The results 
presented in Figure [7] quantify the extent to which the ad¬ 
dition of the Bz component reduces the net torque. The 
mitigating effect of an added vertical field component is 
also evident in the weakening of the density perturbation 
at the MR locations in the B^ + B^ model in comparison 
with the pure-i?,p case in Figures U and [51 

In the left panel of Figure [7] we also present the torque for 
& Bz + Bip field configuration in the case where —>■ 0 
but the condition v'^ = 0 is n ot enforced. As predicted by 
the analysis in Section 12.31 the net torque is somewhat 
larger in this case (implying faster outward migration) 
than in the strict-2D limit. 

3.3.2. Self-Similar Disk 

We carried out similar calculations to those in Sec¬ 
tion |33T] for this case. Our results are presented in Fig¬ 
ure m in which the left and right panels are the analogs 
of Figure El and the right panel of Figure III respectively. 
The semianalytic calculations shown in the left panel of 
Figure [5] demonstrate that, in this case, too, the mag¬ 
netic resonances associated with the presence of a B^^ 
component are shifted away from the corotation radius 
when a Bz component is added. However, in contrast 
with the situation in a uniform, force-free disk, in this 
case rc does not coincide with rp (see Equation (ITT)) and 
Figure |T]), and it is evident from the figure that the reso¬ 
nances are centered approximately on rc rather than on 

rp. 

The density profile shown in the left panel of Figure |8| 
exhibits a noticeable dip at corotation for the pure-i?,^ 
disk model. This feature (which is less than ^1/2 the 
strength of the dips seen at the locations of the MRs 
for this field configuration) could arise from a nonlin¬ 
ear effect since ther e is no linear corotation resonance 
for a pure-R,^ field (|Teraueml 120031 1. Due to its relative 
weakness, this feature does not significantly affect the 
net torque in this example. However, the questions of its 
origin (including why it does not appear in the force-free 
case shown in Figure E]) and of its potential influence in 
other disk models merit further study. 

The right panel of Figure [8] shows that, unlike the situ¬ 


ation considered in Section 13.3.11 the net torque in the 
self-similar case is negative, corresponding to inward mi¬ 
gration. The difference from the behavior of a uniform, 
force-free disk can be understood as follows. We recall, 
first, that in a Keplerian HD disk several factors co mbine 
to ma ke the net torque on the planet negative fe.g.. lWardl 
(HIS), one being the fact that the outer LRs lie slightly 
closer to the planet than the inner resonances of the same 
order. In such a disk, a radial density gradient has lit¬ 
tle effect on the torque since the enhanced weighting of 
a denser inner disk is compensated to a large extent by 
the associated thermal pressure gradient, which shifts r^ 
inward relative to rp and thereby increases the relative 
strength of the outer resonances. In the case of a uniform, 
force-free disk, oc r^, resulting in the outer MR being 
shifted (by magnetic pressure effects) farth er away from 
rc than the inner MR (see Equation (llBbll l. thereby en¬ 
hancing the relative contribution of the inner MR. (There 
is no shift in rc with respect to rp in this case because 
of the absence of thermal and magnetic pressure forces.) 
This has led to the outward migration that we found in 
Section 13.3.11 However, in the present case there is no 
variation in the positions of the inner and outer MRs rela¬ 
tive to rc because is a spatial constant in a self-similar 
disk. At the same time, rc is shifted inward relative to 
rp (as the left panel of Figure El explicitly demonstrates), 
enhancing the relative contribution of the outer MR. As 
a result, the net torque is negative in this case (although 
its magnitude is reduced by ^ 35% when a vertical field 
component of equal amplitude is added to the azimuthal 
field component). This interpretation is also useful for 
under standing the results presented by iFromang et al.l 
(1200^ . In fact, we find that among the different mag¬ 
netized disk models that they explore (listed in their ta¬ 
ble 1), the ones that exhibit outward migration are those 
in which increases with r. (For their remaining mod¬ 
els, whinh_exhibitinwmd migration, fi^p decreases with 
r.) IFromang et al.l ()2005f) discuss the sense of migration 
in their model disks in relation to the radial variation of 
Bp,, but the radial dependence of {3^, is a more pertinent 
quantity in this regard. 

3.4. 3D Simulations 

We now present results from 3D simulations that take 
into account the vertical component of the planet’s grav¬ 
ity (although, unless stated otherwise, we continue to ne¬ 
glect the vertical component of the stellar gravitational 
field). We explore the vertical structure of resonances 
and the effects of wave propagation in the z direction on 
the direction and speed of planet migration. We note 
that our 3D simulations do not have the numerical reso¬ 
lution to separate the MRs from the ARs. However, our 
semianalytic study indicates that the dominant contri¬ 
bution to the density perturbations in this regime comes 
from the MRs. 

3.4.1. Vertical Field in a Uniform Disk 

This case is studied semianalytically in Paper H. Figure|9] 
presents the inidplane density structure from a numerical 
simulation performed using as model parameters a = 0, 
Os = 0, qz = 0, and Pzp = 1- The plot exhibits a dom¬ 
inant one-armed spiral (with an amplitude of ^ 10% of 
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Fig. 7.— Left: Torque on the planet as a function of the azimuthal mode number m, calculated semianalytically for a 2D pure-Bi^ 
field and for a field configuration assuming either fcj = 0 (but with v'^ allowed to be nonzero; 2D) or = 0 (strict 2D). Right: 

Cumulative torque on the planet from the 2D numerical simulations for the same two field configurations shown in the left panel. Both 
calculations employ the same radial scalings (a = 0, as = 0, Qz = 0, qip = 1) and assume /3p = 1 for each field component, although they 
differ in the value of Cp (0.03 in the semianalytic treatment vs. 0.1 in the simulations). The torque units used in the left panel can be 
converted into the units employed in the right panel through a multiplication by (Sprp)/M* fii (2pp/iprp)/M*. 




Fig. 8.— Left: Same as Figure [O] except that the model parameters in this case correspond to a self-similar disk (a = 3/2, as = 1/2, 
Qz = q<p = 5/4). The locations of the corotation radii for the pure-il,^ and Btp + Bz field configurations are also marked (cf. Figure [Q. 
Right: Same as the right panel of Figure[3 but for a self-similar disk model. 


the unperturbed disk density) as well as similar spiral fea¬ 
tures that are, however, much weaker (< 1% of po). The 
stronger spiral resembles the perturbed density struc¬ 
ture of a Keplerian HD disk and can be similarly inter¬ 
preted as a wake formed by the constructive interference 
of acoustic waves of different azimuthal mode numbers 
m th at are launched at the e ffective Lindblad resonances 
(e.g. iQgilvie fc Lubowll2002D . We interpret the fainter 
spiral features as analogous “magnetic wakes” that are 
formed from the constructive interference of SMS and 
Alfven waves launched at the turning points associated 
with the MRs and ARs. In fact, according to Equations 
cm and (fT51) (see also Figure II.l), both ARs and MRs 
are present in this case, and their locations depend on 
m. The locations of the associated turning points are 
likewise m-dependent (Figure II.2). 

As shown by the short-dashes black curve in Figure [TUI 
the net torque in this case is negative — corresponding 
to inward migration — but is reduced by ~ 50% with 


respect to the 3D torque exerted on the planet in a uni¬ 
form, unmagnetized disk (solid black curve). This result 
is qualitatively similar to our finding for this field config¬ 
uration in the 2D limit (Section 13.21 although note that 
we considered a different disk model in that case). One 
can understand this similarity on the basis of our semi¬ 
analytic results (Section II.4.2), which indicate that the 
cumulative 3D torque is comparable to the integrated 
contribution of the kz = 0 vertical mode (i.e., to the 
2D torque) for the pure-Hz, uniform disk model. In Pa¬ 
per II we obtain a rough semianalytic estimate for the net 
torque in this case — a factor of ~ 2 smaller than the 
2D HD torque. We have also estimated semianalytically 
the net torque for an unmagnetized disk in 3D and found 
that is slightly smaller than the 2D HD torqueH Thus, 

^^ ITanaka et al.l ||2002I ~) found that, in the HD case, the 2D torque 
exceeds the 3D torque by a factor of > 2. However, they employed 
different disk parameters, and therefore an exact correspondence 
with our result is not expected. 
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in both 2D and 3D, we estimate that the net torque on 
disks with a pure-Bz field is a factor of ~ 2 smaller than 
the HD torque. This is consistent with our 2D and 3D 
simulation results. 



Fig. 9.— Midplane density in a uniform, 3D disk permeated by 
pure-Bz field with f3zp = 1. The planet is located at r = 1.0 
and ip = TV. The contrast in the figure was adjusted to show, in 
addition to the bright Lindblad wake formed by the constructive 
interference of acoustic (FMS) waves, also the analogous, but much 
fainter, wakes that are formed by SMS and Alfven waves. 


3.4.2. Force-Free Azimuthal Field in a Uniform Disk 

This model (a = 0, Og = 0, q^p = 1) extends th e 2D 
treatment of a pure -i?,„ disk in iTerauemI (|2003ll and 
iFromang et al.1 (j2005[ l to 3D. It incorporates the perturb¬ 
ing vertical gravitational force exerted by the planet and 
the magnetic and thermal-pressure forces arising from 
the disk’s response. 

The behavior of the midplane density for this case is 
shown in Figure fTTl During the first few orbital periods 
of the planet, the density structure appears to evolve 


similarly to that in the 2D simulation (Figure S]). How¬ 
ever, in contrast to the 2D case, the net torque becomes 
negative (dash-dotted black curve in Figure [TOl). The 
reason for this difference is probably similar to the cause 
of the corresponding behavior of the B^ + B ^ field config¬ 
uration, which we discuss in Section [3.4.41 In any case, 
after about 3 to 4 orbits, even before the imprints of the 
resonances have been fully established, the disk devel¬ 
ops an instability roughly at the radial position of the 
the planet. This instability spreads rapidly to the rest of 
the disk (mostly into the inner regions over the duration 
of the simulation), creating turbulence that becomes the 
dominant influence on the subsequent evolution of the 
planet. This evolution is not pursued in this paper as 
it would re quire a signihcant additional computational 
effort (e.g., llJribe et al.l 120111 1. For this reason we have 
truncated the hybrid-field curve in Figure [10] just before 
the onset of the instability. 

A more detailed characterization of the observed insta¬ 
bility is as follows. Initially, the planet creates small per¬ 
turbations around its orbit in all the physical quantities 
(density, velocity, and magnetic field). The vertical field 
component starts growing after about 1.5 orbits, and the 
radial field component follows suit after ~ 2.5 orbits (see 
left panel of Figure nil). The two poloidal field com¬ 
ponents then grow exponentially together, increasing in 
energy by ^ 4 — 5 orders of magnitude until they satu¬ 
rate after ^ 6 orbits. The evolution of the associated (3 
parameters is shown in the right panel of Figure [lH The 
figure also shows that the energy associated with the az¬ 
imuthal field, which is dominated by the ordered (equi¬ 
librium) component of B^,^ starts decreasing when the 
instability is triggered, with the decline becoming more 
pronounced after the instability saturates. We have ver¬ 
ified that this decline is associated with the radial mass 
inflow that is induced by the onset of MRI turbulence: 
the large-scale magnetic field is then simply advected in¬ 
ward by the inflowing mass and leaves the computational 
domain through the inner radial boundary. We did not 
attempt to address this issue in our simulations because, 
on the one hand, it does not affect our conclusions about 
the development of the instability (whose growth occurs 
on the dynamical timescale, which is much shorter than 
the radial inflow time), and, on the other hand, we do 
not attempt to model the longer-term evolution of the 
disk accurately. 

The evolution of the resulting turbulent stresses in the 
disk is shown in Figure 1131 They start growing expo¬ 
nentially after ^ 4 orbits, and saturate after ~ 6 orbits. 
It is seen that the Maxwell and Reynolds stresses reach 
comparable magnitudes (with the Maxwell stress being 
larger overall). The total stress (normalized by the ini¬ 
tial thermal pressure at and presented as an effective 
viscosity parameter) grows to a ~ 0.lE3 

To check whether the instability might be stabilized by 
vertical density stratification, we ran a simulation with 
the same numerical configuration except that we also in¬ 
cluded the vertical component of the stellar gravity. We 

We consi der only t he rip component of the stress and use the 
expressions in IHawlevI 1120001 ) for the vertically and azimuthally 
averaged Reynolds and Maxewell contributions at any given radial 
position. 
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Fig. 10.— Cumulative torque for all of the reported 3D simulations of an unstratified disk. 


found that the instability still developed in the result¬ 
ing stratified disk, and that the magnitude of the torque 
exerted on the planet before the development of the in¬ 
stability was not significantly different from the value 
obtained for the unstratified disk. The difference be¬ 
tween these two cases is illustrated in Figure [Ml which 
shows the poloidal (r — z) density structure of the MRs 
that develop inside and outside the planet’s orbit before 
the instability is triggered. The left panel of this figure 
demonstrates that the radial position of each of the MRs 
is constant with height for the uniform disk, z, whereas 
the right panel shows the fanning-out of the resonance in 
the stratified disk. This behavior is a dire ct consequence 
of the dependence of tmr (Equation (IMbp i on and the 
fact that I3ip oc p decreases with z in a stratified disk. The 
observed density structure also reflects the weakening of 
the perturbing gravitational potential with distance from 
the planet, which causes the resonance amplitude to de¬ 
cline with height. 

We discuss the nature of this instability in Section 01 
where we argue that it is essentially a form of the MRI. 
The important implication of this result to planet migra¬ 
tion is that the conclusion s drawn for this field configura¬ 
tion f rom 2D studie s fe.g.. lTerauemll20r)^lFromanef et al.l 
120051 Section 13. 3. IL and in particular the prediction of 
a systematic outward migration for large enough values 
of are not applicable to real disks, which are inher¬ 
ently 3D. A planet embedded in a disk that is perme¬ 
ated by a pure-Ri^ held would evolve, instead, according 
to the predictions of MRI -turbulent disk models (e.g., 
iNelson fc PapaloizorjDOO^ . In particular, a planet with 
the mass employed in our simulations (~ would 

experience the torque responsible for Type I migration 
in a laminar disk as well as a measurable random torque 


component induced by the turbulence, likely resulting 
in a nonmonotonic, but overall inward-directed, drift. 
The relative impact of the stochastic component is pre- 
dicted to increas es with decreasing planet mass (e.g., 
lUribe et al.ll2011l) . 

3.4.3. Disk with B^p oc z 

We now consider the same disk model as in Section l3.4.2l 
except that, instead of having a force-free azimuthal held 
that is uniform in z {Bq^ = we add a depen¬ 

dence on the vertical coordinate: Bq^ = B^p^^z/ZQ)r~^, 
where zq = 0-13 is the upper domain boundary (see Sec¬ 
tion |3T|^_^e value of Bp,p is the same as that used in 
Section 13.4.21 and corresponds to /3c^p = 1 if one uses the 
midplane density in evaluating this parameter. (The ac¬ 
tual value of /3ip aX z = zq is, however, < 1 since the den¬ 
sity there is lower than at the midplane.) The adopted 
vertical structure of By, is intended to approximate the 
azimuthal held in a wind- driving disk that is thr eaded 
by open held lines (e.g., iWardle fc Konigll 1199311 . and 
we consider it here (as well as semianalytically in Sec¬ 
tion 1.4.3) in order to isolate the effects of the dBy,/dz 
term on the disk behavior. The main expected impact 
would come from the magnetic torque term in the an¬ 
gular momentum equation {{B^/pi)dBy,/dz) — however, 
this term (which gives rise to radial motion that, in turn, 
generates a B^ component in a disk threaded by open 
held lines) cannot be included in an ideal-MHD simu¬ 
lations that aims to approach a steady state (see Sec¬ 
tion [TJQ 

The z-gradient term also appears in the vertical momen- 

We are, however, able to infer the effect of the magnetic torque 
term on the linearized equations; see Paper II. 
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Fig. 11.— Snapshots of the midplane density (labeled by the number of planet orbits) in a uniform, 3D disk with a force-free field. 
The planet is located at r = 1.0 and ip = tt. The development of an instability after ~ 3 — 4 orbits is clearly visible. 


turn equation {—{d/dz)[B‘^/{2fi)]), where it represents 
magnetic pressure squeezing that has the same effect as 
the tidal squeezing by the stellar gravity. One can there¬ 
fore anticipate that the poloidal structure of the reso¬ 
nances in this case will be similar to that of the stratified 
disk model considered in Section [3.4.21 This is confirmed 
by the results shown in Figure [121 which resemble those 
in the right panel of Figure HT] We initialized the sim¬ 
ulations presented in Figure [15] using the approximate 
vertical-equilibrium expression derived in Appendix m 
the density adjusts to a fully consistent steady state af¬ 
ter a few orbits. We show results for a model in which 
just the magnetic squeezing effect is included as well as 
for a disk in which both magnetic and tidal squeezing op¬ 
erate. According to Equation (IB2I) . the ratio of the mag¬ 
netic and tidal squeezing terms is ^ {hp/ zq)"^, which 
equals 0.59 for our chosen parameter values —- indicat¬ 
ing that the magnetic squeezing associated with the By,- 
gradient term has only a small effect. This is consistent 
with more realistic mod els of wind-driving disks (e.g., 
IKonigl &: Salmer'onll2011l) . which indicate that magnetic 
squeezing is important in well-coupled disks but is in¬ 
duced by the radial, rather than the azimuthal, field 
component. (As we noted above, we cannot include a 
Br field in an ideal-MHD simulation of this problem.) 

The fact that the amplitude of the azimuthal field compo¬ 
nent remains small in the vicinity of the planet weakens 
the effect of the resonances, which only develop above 


the midplane in this case. Consequently, the torque on 
the planet in this model is effectively the same as in an 
unmagnetized disk (see dash-triple-dotted black curve in 
Figure [ini where an unstratified disk is considered). This 
result is consistent with the inference from our semiana- 
lytic study of this field configuration (with = 0 instead 
of 1) in Paper II. Another consequence of the small ampli¬ 
tude of the equilibrium magnetic field near the midplane 
is that the perturbed vertical magnetic field is evidently 
too weak (for the given level of numerical diffusivity) to 
trigger the MRI, so that, in contrast with the behavior 
of the uniform-disk (see Sections 13.4.21 and [4]) , no 
instability develops in this case. 


3.4.4. Force-Free By + Field in a Uniform Disk 

We now combine the vertical and azimuthal field com¬ 
ponents discussed in Sections 13.4.11 and 13.4.21 respec¬ 
tively, and consider a hybrid B^ + By disk model with 
a = as = 0, qz = 0, Qy = 1, and Pzp = (3yp = 1. For 
this field geometry, both ARs and MRs are excited — 
their locations given by Equations m and m , respec¬ 
tively — and the associated waves (SMS and Alfven) 
propagate along the direction of the total field. Since 
the resonances and turning points depend on m (see the 
bottom panel of Eigure[3|), the density perturbations will 
again form magnetic wakes, analogous to those found in 
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Fig. 12.— Time evolution of the magnetic energy associated with each field component (Left) and of the corresponding /3 parameters 
(Right), spatially averaged over the full computational domain, for a 3D, uniform disk with a force-free field. 
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Fig. 13. — Time evolution of the Reynolds, Maxwell and total 
stresses, each normalized by the initial thermal pressure at the 
planet’s radial di stan ce from the star (1 AU), for the model con¬ 
sidered in Figure ^2] The stresses are averaged over the azimuthal 
and vertical directions and are evaluated at a disk radius of 0.8 AU. 

the pure-B^ case (Figure 

The total torque on the planet in this case is given by the 
black long-dashes line in Figure [TOl Just as in the case 
of each of the separate components (with the pure-i?,^ 
disk considered prior to the onset of the instability), the 
magnitude of the torque for the hybrid field configuration 
remains smaller than the HD value. However, in contrast 
with the 2D case (see Figure ED, it is larger (although 
not by much) than the magnitude of the pure-i?,^ torque. 
Since the magnitudes of the 3D torque contributions from 
the inner and outer disk regions are both larger in the 
pure-H,^ model, the dominance in 3D of the net torque for 
the hybrid-field configuration is a consequence of a shift 
in the balance between these two contributions on going 
from 2D to 3D. The origin of this shift may be related 


In practice, the ARs are not visible or even resolved in the 
numerical simulations since the MR-to-AR amplitude ratio is typ¬ 
ically ~ 10"* for the pure-S,^ case and ~ 10® for the hybrid-field 
case. 


to the apparent asymmetry in the radial extent of the 
innermost evanescent region between the inner and outer 
disk regions in the 3D hybrid-field case, particularly for 
low values of the azimuthal mode number m (see lower 
panel of Figure [31). Note that our turning-point analysis 
in Section r2.3l also indicates that the size of the innermost 
evanescent region in the m — r plane is smaller, in both 
the inner and outer disk regions, for the 3D pure-H,^ 
model (upper panel of Figure|3]), which is consistent with 
our finding that the contribution to the torque from each 
of these regions is larger in this case than in the hybrid- 
field model. 

A notable difference between the cumulative 3D torque 
for this field configuration and the 2D torque presented 
in Figure|7]is that in 3D the total torque on the planet is 
negative, corresponding to inward migration, whereas in 
2D it is positive (outward migration). A similar differ¬ 
ence was already seen in the behavior of the pure-H,^ field 
(again considering the situation in the 3D case before the 
onset of an instability). To understand the origin of this 
difference, we examine the contributions to the torque 
from different regions of the disk based on our linear 
analysis results. Tabled] shows the torque exerted on the 
planet — for particular azimuthal and (in 3D) vertical 
mode numbers — by regions outside the LRs and by re¬ 
gions surrounding the magnetic and Alfven resonances 
for both 2D and 3D realizations of this model. 

In the case of a pure-5,^ field in 2D, iTerauenJ l|2003[) 
found that the dominant contribution to the torque 
comes from the regions around the MRs even though the 
torque from the point-like contribution at the resonances 
can sometimes be of opposite sign. We have confirmed 
this result and found that it applies also to the other 
cases considered in Table [l] The point-like contribution 
from the MR can be opposite of that from the surround¬ 
ing region for to < 40 (although not for higher values 
of to), and the overall MR contribution dominates the 
total torque. In 2D, the positive contribution from the 
inner MR region has a larger magnitudes than the nega¬ 
tive contribution from the outer MR region for both the 
purc-i?,p and the + B^p field configurations, resulting 
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Fig. 14.— Polo ida l de nsity structure for two 3D disk models with a force-free azimuthal field before the development of the instability 
shown in Figures Il2l and [13I The left panel represents a uniform disk, whereas the right panel shows a disk that is vertically stratified by 
the action of the stellar tidal gravity. The planet is located at {r,z) = (1,0). 



Fig. 15.— Same as Figure fT^ but for an azi muthal field that increases linearly with height from zero at the midplane to its final value 
(which equals its constant value in Figure fT^ at the top of the computational domain. The vertical gradient of contributes to the 
vertical stratification of the disk. The disk shown in the right panel includes the effects of both the field gradient and the stellar gravity 
and is therefore more strongly stratified than the model presented in the left panel, in which the stellar contribution is neglected. This is 
manifested in the larger opening angle of the MR “perturbation cone” in the right panel. 


in a positive net torque on the planet 

For the 3D hybrid case, this trend is reversed: The outer 
MR region dominates the inner one, resulting in a nega¬ 
tive net torque. This behavior can be understood from a 
consideration of the MR locations, using Equations m 
and (US for the 2D and 3D cases, respectively. Whereas 

The Tio torque value given in Table^for the 2D pnre-Bp disk 
model (73 in the adopted units) c an be directly com pared with the 
corresponding entry in table 1 of ITeroueml 112003) . which is 66 in 
the same units (taking account of the opposite sign convention 
employed in that table). The numerical difference between these 
two values reflects the difference in the adopted radial range of 
integrati on (0.3 — 1.7 rp in our semianalytic calculations vs. 0.5 — 
1.5 rp in ITeroueml 12003 : we reproduce the value in table 1 of the 
latter reference when we integrate over the narrower range). 


the outer MR in the force-free case always lies farther 
away from the planet than the inner MR, in 2D the dis¬ 
tances of the inner and outer MRs from rp are compa¬ 
rable (their difference is ^ 1% of their mean). When 
kz ^ 0, the asymmetry between the distance of the outer 
MR and inner MR is increased (the difference is ^ 18% 
of the mean for kzh = 1.56). This leads to the point-like 
contribution at the outer MR being significantly weaker 
relative to its surrounding region in the 3D case than in 
2D. Hence, in 3D, the outer point-like contribution sub¬ 
tracts less from the torque arising from the surrounding 
region, resulting in that region’s contribution dominat¬ 
ing the net torque. As is seen from Table [U the effect 
of the MR regions in 3D is reinforced by that of the AR 
regions, which also contribute a negative torque from the 
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outer side that outweighs the positive contribution of the 
inner side. As a check on these inferences, we decom¬ 
posed the numerically derived density for the hybrid-field 
case into Fourier modes in the vertical direction, using 
kz = ^imz/Lz (where Lz is the vertical domain size and 
Uz is a nonnegative integer). We determined that the 
sum over the > 0 modes dominates the Uz = 0 am¬ 
plitude in the region around the outer MR. This is con¬ 
sistent with the linear-analysis arguments given above 
and supports the conclusion that the kz ^ 0 modes are 
responsible for the reversal of the outward migration in¬ 
dicated in the 2D limit. 

Unlike the pure-B,^ disk, our hybrid field model (charac¬ 
terized by equal amplitudes of the azimuthal and vertical 
equilibrium field components) remains stable in 3D. We 
discuss the reason for the difference between these two 
cases in Section IH 

4. DISCUSSION 

Our ultimate goal is to model planet migration in a 
generic wind-driving disk where the magnetic field is ev¬ 
erywhere well-coupled to the matter. The midplane mag¬ 
netic field configuration of such a disk is characterized 
by Bz ^ 0, = Br = 0, dB^/dz < 0, dBr/dz > 0, 

whereas near the surface the amplitudes of all three field 
components are comparable (with \Br\-, however, typ¬ 
ically exceeding |B,^|). A defining property of such a 
disk is that the magnetic torque term oc BzdB^p/dz ac¬ 
counts for a significant fraction — and possibly all — of 
the angular momentum transport that enables local gas 
to flow toward (and ultimately accrete onto) the cen¬ 
tral star. This transport is vertical, and in the case of 
a wind-driving disk the angular momentum is deposited 
in an outflow driven from the disk surfaces. As we have, 
however, pointed out, this model cannot be studied nu¬ 
merically with an ideal-MHD code, and a linearization 
treatment is similarly hampered by the need to incorpo¬ 
rate magnetic diffusivity in order to establish an equi¬ 
librium field configuration in the differentially rotating 
disk. For this reason, we have concentrated in this paper 
on model problems that can be investigated within the 
framework of ideal MHD. While these problems do not 
directly mimic the expected situation in a wind-driving 
disk, they can increase our understanding of the possible 
effects of a large-scale, ordered magnetic field on Type I 
planet migration. It is also conceivable that some of the 
field configurations that we explore might actually apply 
to real systems in certain circumstances. 

The models that we considered — pure-Bz-, puie-B,^, and 
a hybrid Bz + configuration, in either a force-free or 
a self-similar disk — were examined in the 2D {kz —>■ 0) 
and strict-2D (u(, —>■ 0) limits (with the latter case cor¬ 
responding to the 2D implementation in numerical sim¬ 
ulations) as well as in 3D using both semianalytic and 
numerical techniques. We now discuss the main insights 
provided by this study, which include the effect of a ver¬ 
tical field component, the difference between 2D and 3D 
models, and the onset of a disk instability in 3D. 

As was pointed out by iMuto et all (1200811 and confirmed 
in our work (see Sections [Ql and 13.21 a purely verti¬ 
cal field in 2D simply acts to increase the effective sound 
speed in the disk, which pushes the turning points farther 


away from the planet and enhances the sharp, pressure- 
induced cutoff in the torque, thereby reducing the mag¬ 
nitude of the net torque acting on the planet in com¬ 
parison with the unmagnetized disk case. We found that 
this is also the main effect of adding a Bz component to a 
disk in which an azimuthal field is present (Section [TS]): 
the magnitude of the net torque is reduced but its sign 
(positive or negative) remains unchanged. Thus, in a 
force-free disk {B,^{r) oc r“^) for which outward migra¬ 
tion is predicted, the addition of a uniform vertical field 
Bz = \B,p{rp)\ reduces the cu mulativ e torque by ~ 35% 
but keeps it positive (Section 13.3.11) . whereas in a self¬ 
similar disk {B,^ oc r“®(^), for which inward migration is 
inferred, the addition of an equal-amplitude vertical field 
component reduced the magnitude of the net torque by 
^ 50% but keeps it negative (Section 13.3.21) . This be¬ 
havior can be understood from the fact that a vertical 
field component does not excite any new resonances in a 
2D disk: its only impact is through the increase in the 
effective sound speed. The consequences of the increase 
in Ceff depend, however, on the nature of the disk model: 
For example, in a 2D model the addition of a vertical 
component to an azimuthal field shifts the MRs toward 
the planet, whereas in a strict-2D model such an addition 
pushes them away (Section [521) • 

In contrast to its 2D behavior, a Bz field in 3D gives rise 
to resonances — both MRs and ARs (with the former 
being stronger than the latter), resulting in the vertical 
propagation of SMS and Alfven waves and the establish¬ 
ment of weak magnetic wakes (although the Lindblad 
wake, associated with the propagation of FMS waves that 
are launched at the effective Lindblad resonances, still 
dominates the density perturbation; see Figure E]). The 
opening up of additional wave propagation channels can 
be expected t o incr ease the torque on the planet, yet we 
find (Section 13.4.11) that the net torque is, in fact, re¬ 
duced (by ^ 50% relative to the HD case for a uniform 
disk with = 1). Based on the results of a linear analy¬ 
sis (Section II.4.2), we attribute this outcome to the fact 
that in this case the 2D mode dominates the contribu¬ 
tions of the kz ^ 0 vertical modes to the net torque. 

As we have seen, a pure-H^, disk that harbors a planet 
becomes rapidly unstable in 3D, but the addition of a ver¬ 
tical field component of comparable strength stabilizes it. 
The hybrid field configuration also develops (dominant) 
MRs and (subdominant) ARs, with the SMS and Alfven 
waves now propagating along the total field. Assuming a 
force-free field in a uniform disk, we find (Section I3.4.4p 
that the net torque in this case is negative — the con¬ 
verse of the situation in the strict-2D limit — which we 
argue is a consequence of a fc^-dependent asymmetry be¬ 
tween the locations of the inner and outer MRs. It is 
noteworthy that the 3D torque for a force-free pure-H,^ 
field (measured before the disk becomes unstable) is also 
negative, and that this is again the converse of the situ¬ 
ation for the corresponding 2D case. This suggests that, 
in 3D, the role of the Bz component in the hybrid field 
configuration is secondary to that of the B,^ component, 
just as it evidently is in 2D. It is also worth noting that 
the magnitude of the 3D hybrid-field torque, while re¬ 
maining smaller than in an unmagnetized disk, exceeds 
the magnitude of the (pre-instability) pure-i?,^ torque. 
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TABLE 1 

Torques from the Different Regions of a Force-Free, Uniform Disk in 2D and 3D 



2D (for 

m = 10) 




B Tu 

TmRi 

TmRo 

Tlo 

^total 


Bz + 39 

257 

-233 

-45 

18 


Bp 76 

507 

-416 

-95 

72 


3D (for 

m = 10 and kzh ■■ 

= 1.56) 



B Tii TARi 

TmRz 

TmRo 

Taro 

Tlo 

^total 

Bz + Bp 1 114 

316 

-422 

-134 

-1 

-126 


Note. — In 2D, the torque regions are defined as follows: is the torque from the region extending from the outer edge of the disk 

to Rl+, and Tmro is the region surrounding the MR from Rm+ to In 3D, and Tmro enclose T/^r^ , which we define as the 
region surrounding the outer AR going from Ra+ to Ra-- See Figures|2]and |3]for the turning-point labels. The inner regions (subscript 
‘i’ instead of ‘o’), are defined analogously. 


which contrasts with the situation in the 2D hybrid-held 
models (Section 13.31) . This difference might also be the 
result of a subtle effect such as turning-point asymmetry, 
and it provides another illustration of the fact that the 
details of the behavior in 3D are in general more complex 
(and less predictable in advance) than in 2D. 

A potential caveat to our results is that the hnite verti¬ 
cal domain of the simulations constrains the modes that 
can be resolved numerically; only modes with wavenum¬ 
bers kz > kz^mi-a = 27r/Lj; = 24.17 (for = 0.26) 
can play a role in the simulation. This corresponds to 
kz,minh = 2.42. In the case of unstratihed disk models, 
Lz is the only lengthscale that influences the relative im¬ 
portance of the vertical modes, which in principle can 
lead to an unphysical dependence of the simulation re¬ 
sults on the domain length. We have verified that this 
is not a significant effect in our simulations by running 
vertically stratified disk models for the uniform-and 
B^p{z) oc z field configurations and determining that the 
calculated torques are very similar to those obtained for 
the unstratified disk models. We also doubled the ver¬ 
tical domain size for the simulated uniform hybrid-held 
model and confirmed that the derived torque remained 
essentially the same. 

The pure-i?,p field configuration considered in Sec¬ 
tion 13.4.21 exhibits an instability that develops at the 
position of the planet and, within a few orbital periods, 
results in strong turbulence that continues to spread into 
the surrounding disk. The characteristics of this insta¬ 
bility — including the growth of the turbulent Reynolds 
and Maxwell stresses on a timescale on the order of 
with the stresses remaining comparable (yet dominated 
by the magnetic contribution) and saturating to a total- 
stress-to-gas-pressure-ratio of a > 10“^, as well as the in¬ 
sensitivity to background stratification — have the hall¬ 
marks of a “standard” MRI , which arises in a Keple- 
rian disk that is thread ed by a weak vertical field (e.g., 
iBalbus h Hawle^ll998[) . This interpretation may at first 
seem puzzling since a similar instability does not de¬ 
velop for either the pure-R^ or the hybrid {B^ + Bz) 
field configuration. The absence of this instability in 
the latter two models can, however, be understood from 
the existence of a minimum vertical wavelength for the 
development of a standard MRI (arising from the need 
to overcome the stabilizing effect of the magnetic ten¬ 
sion force) and from the requirement that this wave¬ 


length remain less than the disk thickness. This condi¬ 
tion can be expressed as a lower limit on the parameter 
Pz'- Pz,ram = (tt/a/S)(/ i/ro) for instability. For /ip = 0.1 
and Zq = 0.13, the values used in our simulations, we 
obtain /3z,min = 1-4. Since our pure-Rz model employs 
Pzp = 1, and so does also our hybrid model, both are 
stable to the standard MRI. 

A pure-R;^ disk is also subject to a weak-field shear¬ 
ing instability, although, unlike the standard MRI, 
it is triggered by the excitation of non axisymmetric, 
rathe r than axisymmetric , modes (e.g., iHawlev et al.l 
1199,^ iRiidiger et al.l 120071 ). In the limit of large az¬ 
imuthal wave numbers (m > h ~^), one can write down 
the instability condition (e.g ., iQgilvie fc Pringld 119961 
iTerauem fc Panaloizod [T9M ) and again convert it into 
a lower bound on the relevant P parameter: P^p > 
Pip,rain = mh/\/i for instability. The numerical value of 
/3c,5,min should be comparable to that of /3z,min, and one 
can therefore expect that our pure-R,^ disk model, with 
Pp-p = Ij would also be stable to the MRI (the azimuthal 
MRI (AMRI) in this instance). In fact, even in the case of 
instability, the growth rate of the AMRI would be much 
slower than that of the standard MRI (i.e., -C D“^) un¬ 
less kzh were ^ 1, in which case the presence of even a 
tiny poloidal field co mponent could transform the nature 
of the problem (see IBalbus fc HawlevI 1199^ . It seems 
that this is precisely what has transpired in our 3D pure- 
Bp simulation; This configuration evidently remains sta¬ 
ble to the AMRI, but the growth of the perturbed verti¬ 
cal field component (of initial amplitude < 10~^\Bop\), 
induced by the gravitational perturbation of the planet, 
triggers the standard MRI and turns the disk turbulent 
over a few orbital period s. We note in this connection 
that IHawlev et al.1 (jl995f ) previously demonstrated that 
the MRI behavior of a hybrid B^ -|- Bz field configuration 
is roughly the same as that of the sum of the two indi¬ 
vidual models; in particular, they found that the evolu¬ 
tion of a standard MRI triggered by a weak vertical field 
is unaffected by the presence of a comparatively strong 
(AMRI-stable) azimuthal field in the disk. 

The field configurations that we considered are po¬ 
tentially subject also to other types of instabili- 
ties. For example, the h ybrid case was shown by 
iHollerbach fc Rudigeil (|2005f ) to be susceptible to the de¬ 
velopment of a nonstationary flow pattern, consisting of 
elongated Taylor vortices that drift in the vertical di- 
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rection. This instability can be attributed to the hand¬ 
edness of this field configuration, which is only present 
when both the and B^p components are nonzero. In 
our simulation of the hybrid field model we also observe 
the development of a nonsteady flow pattern that breaks 
the ±z symmetry, but a longer run that covers a larger 
domain is needed to study this behavior in full detail. 

5. CONCLUSION 

We have studied Type I migration in protoplanetary 
disks that are threaded by a large-scale, ordered magnetic 
field. Among the possible origins of such a field in real 
disks are interstellar field lines that have been carried in 
by the accretion flow and a dynamo-generated field from 
the central star that has diffused into the disk gas. In ei¬ 
ther of these two cases, the field can be modeled as being 
open and possessing an even symmetry (Br = B^ = 0, 
i ?2 0 at the midplane). In general, the drag exerted 

by the inflowing gas and the differential rotation of the 
disk give rise to vertical field gradients — {dB^/dz and 
dBp/dz — that, in conjunction with the vertical field 
component B^-, produce a radial magnetic tension force 
and a braking torque, respectively. The torque term can 
lead to significant vertical transport of angular momen¬ 
tum, which, in turn, may be deposited in a powerful 
wind that is driven centrifugally from the disk surfaces. 
Strong outflows of this type are a ubiquitous feature of 
protoplanetary systems, and it is of interest to explore 
what effect the magnetic angular momentum transport 
associated with them might have on planet migration in 
the wind launching region. However, the general prob¬ 
lem cannot be properly studied within the framework 
of ideal MHD because the strong winding up of the field 
lines in a Keplerian disk leads, in the absence of magnetic 
diffusivity, to an ever growing torque that results in an 
unchecked growth of the inflow velocity. A further com¬ 
plication is associated with the fact that a wind-driving 
disk cannot be treated as a spatially localized problem. 
For these reasons we have chosen to concentrate in this 
paper on simplified field configurations that can be han¬ 
dled in a restricted domain and using ideal MHD, in the 
hope that they could provide useful insights into the more 
general problem. Some of these field geometries might 
even be relevant to real systems, and, in any case, their 
treatment enables us to make contact with (and gener¬ 
alize) previous studies in the literature. 

The field geometries that we considered are: pure-H,^ 
(previously stu died semianalyt icall y and numerically in 
the 2D limit bv lTerauemll20Ci3l and lFromang et ^120051 
respectively), pure-i?^ (previously considered semiana¬ 
lytically and numerically, but exclusively in the shearing- 
sheet appro ximation and emp loying additional simplifi¬ 
cations, bv iMuto et al.l 1^0811 . and a hybrid {B^ + B^) 
field configuration, where the fiducial (3 parameter value 
for each of the field components was, in all instances, I. 
In the case of the pure-H,^ field, we studied both a ver¬ 
tically uniform configuration and a model in which the 
field grows linearly with height from zero at the mid¬ 
plane. We examined uniform, force-free disks, in which 
neither a thermal nor a magnetic force plays any role in 
the equilibrium state, as well as disks in which the radial 
dependence of the physical variables corresponds to that 
of the radially self-similar model of wind-driving disks. 


Our study comprised a linear perturbation analysis and 
numerical simulations in both 2D and 3D, where we dis¬ 
tinguished between the analytic 2D limit of setting the 
vertical wavenumber kz equal to zero and the stricter 
numerical limit of letting the vertical component of the 
perturbed velocity vanish. Our linearization analysis has 
yielded the relevant resonances for the different field con¬ 
figurations and disk geometries and the associated wave 
propagation regions. We identified magnetic and Alfven 
resonances and obtained their locations with respect to 
the corotation radius, and we presented algebraic expres¬ 
sions for these quantities (derived under the assumption 
of z invariance). 

For any given model, we derived semianalytically the 
turning-point locations that delineate the propagation 
regions of the respective waves (slow-magnetosonic and 
Alfven) as well as the locations of the magnetically 
modified effective Lindblad resonances from which fast- 
magnetosonic waves are launched. We then evaluated 
the torques contributed by different regions in the disk 
and attempted to estimate the magnitude of the cumu¬ 
lative torque for that case. This information was used 
to analyze and complement the density maps and the 
torque running-time averages obtained in the numerical 
simulations. The main findings of this study can be sum¬ 
marized as follows: 

• The inference from 2D studies of pure-H,^ disk 
models that a sufficiently strong outward decrease 
in the magnetic field amplitude would lead to out¬ 
ward planet migration does not apply to the 3D 
models that describe real systems. In 3D, planet 
migration is invariably inward for each of the field 
configurations and disk models that we considered. 

• In the presence of a planet, a disk threaded by a 
uniform azimuthal field that is characterized by 
a field-strength parameter large enough to 
forestall the development of the azimuthal MRI 
(AMRI) becomes unstable at the planet’s location 
to the standard MRI (which is associated with the 
presence of a weak vertical magnetic field) and 
turns turbulent within a few orbital periods. This 
outcome is found irrespective of whether the effect 
of the vertical component of the stellar gravity is 
taken into account. The instability is triggered by 
the growth of the perturbed vertical magnetic field 
component that is induced by the disk’s interaction 
with the planet. This instability did not develop in 
either the pure-i32 or the hybrid field model since 
our chosen value of I3z was large enough to forestall 
the standard MRI in both of these cases. 

• Even though the presence of a large-scale field can¬ 
not (in 3D) reverse the inward migration predicted 
by the underlying unmagnetized disk model, it gen¬ 
erally slows it down — by a factor of ~ 2 for a 
pure-H^ configuration and a somewhat lower value 
for the hybrid-field case (using Pzp = Ppp = I). 

In the companion paper we carry out a linear perturba¬ 
tion analysis of a disk model that incorporates the ver¬ 
tical magnetic angular transport term (oc BzdB^jdz) in 
an attempt to get a more direct handle on the planet 
migration properties of wind-driving disks. We plan to 
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eventually also perform a 3D nonideal-MHD simulation 
to gain even stronger insights into this problem. 
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APPENDIX 

A. SOFTENING PARAMETER IN 2D SIMULATIONS 

It is well known that the derived magnitude of the torque in 2D simulations depends on the choice of softening 
parameter e in Equation (EOl). For disk-planet interactions in unmagnetized disks, it has been argued that the 
softening parameter should be about 0. 6 h in order for the torque values to match the ones obtained in 3D simulations 
(jPaardekooper fc Panaloizoul l2009a] lbl: iPaardekooner et al.l l201fl[ l . In this work we do not use this scaling for the 
following reasons. First, we do not know whether this particular choice of e is also applicable to magnetized disks. 
Second, the magnetic resonances that play a key role in 2D magnetized disk models generally lie closer to the planet 
than the Lindblad resonances that dominate in the unmagnetized case, which adds the constraint that the distance 
modeled by this parameter should be less than the radial distance from the planet to the MRs. Third, we do not 
attempt to match the 2D and 3D torque values since we are actually interested in the difference between the 2D and 
3D disk models: In the magnetized case, new resonances appear in 3D, along with new regions and directions of wave 
propagation, making the 3D problem fundamentally different from the 2D one. Now, the 2D Lindblad and corotation 
contri butions, which are the relevant torq ue term in the HD case, were shown to converge to a given value when e is 
small (|Paardekooper fc Papaloizoull2009bll . so we take this value to represent the 2D HD torque in the present study. 
Motivated by these considerations, we fix the magnitude of e at 0.03, which is small in comparison with h but still 
large enough to prevent the occurrence of numerical divergences in the vicinity of the planet. 


B. INITIALIZING THE DISK DENSITY PROFILE IN THE PRESENCE OF VERTICAL STELLAR GRAVITY 


We derive an approximate solution for the vertical density structure of a disk that is subject to compression by the tidal 
gravitational field of the central star as well as by the magnetic pressure gradient associated with the azimuthal field 
component. This solution is used to initialize the simulations presented in Figure [15] in order to avoid the development 
of artificial currents. For a vertically isothermal disk, the 2 component of the momentm equation o reads 

2dp GM B^dB^ , 

^ “ ^(r2-Fz2)3/2^ ^TT dz ' ^ ^ 

Under our adopted assumption, B^ oc z (see Section [3.4.3l) . so we can write B^ = {dB^p/dz)z, with dB^jdz = B^p^jz^ 
a constant. Assuming a geometrically thin disk (z r) and focusing on the radial location of the planet, the 
approximate solution of Equation (ED is 


Po(^) 


Pp exp ■ 



(B2) 


where is the (normalized) tidal scale height fSection lXT]) and /3^ is evaluated at the top of the disk (where B^p = H^p) 
but using the midplane density Pp (an approximation justified by the fact that the density does not change much over 
either hp or zq for the adopted values of these parameters ). 
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